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ABSTRACT 

The  extension  of  time-marching  procedures  to  low  Mach  number  and  low 
Reynolds  number  conditions  is  considered.  It  is  shown  that  the  disparate 
speeds  of  the  acoustic  and  particle  waves  prevents  convergence  at  high 
Reynolds  numbers  while  the  requirement  that  both  the  Courant  and  the  von 
Neumann  numbers  be  of  order  one  prevents  convergence  in  very  viscous 
flows.  A  perturbation  expansion  is  used  to  introduce  pseudo-acoustic 
waves  that  propagate  at  speeds  similar  to  the  particle  speed  at  high 
Reynolds  numbers  and  that  allows  both  the  inviscid  and  viscous  time  step 
parameters  to  be  of  order  one  at  low  Reynolds  numbers.  The  resulting 
algorithm  is  shown  to  give  convergence  rates  that  are  independent  of 
either  Mach  number  or  Reynolds  number  over  a  range  of  five  orders  of 
magnitude  in  both  parameters.  Results  are  shown  for  strong  heat  addition 
in  low  speed  flow  encompassing  this  broad  range  of  variables. 

INTRODUCTION 

Time-dependent  algorithms  are  nearly  the  exclusive  choice  for  the 
computation  of  compressible  flows.  They  have  been  highly  developed  to 
apply  to  high  speed  flows  in  general,  and  to  deal  with  the  shock  waves 
that  frequently  appear  under  such  conditions,  in  particular.  Both 
explicit  and  implicit  procedures  have  been  used  extensively  in  transonic, 
supersonic,  and  hypersonic  regimes.  An  important  advantage  of  these 
algorithms  is  that  they  provide  accurate  predictions  in  both  inviscid 
flows  and  in  the  practically  important  regime  of  high  Reynolds  number 
viscous  flows.  One  reason  for  this  flexibility  is  that  they  allow  the 
convective  terms  to  be  central  differenced  at  all  Reynolds  numbers.  In 
cases  where  central  differences  are  not  desired,  they  provide  a  physical 


basis  for  defining  various  upwind  differencing  schemes  also.  In  general, 
upwind  differencing  becomes  more  desirable  as  the  Mach  number  increases. 

An  important  attribute  of  any  computational  algorithm  is  that  it  be 
robust  over  a  wide  range  of  flow  conditions.  In  this  regard,  a  major 
drawback  of  time-dependent  procedures  is  their  well-known  inefficiency  at 
low  subsonic  speeds.  This  characteristic  can  lead  to  difficulties  in 
computations  of  transonic  flowfields  that  contain  embedded  low  speed 
regions  such  as  near  a  stagnation  point  or  in  the  boundary  layer.  In 
addition,  this  limitation  makes  the  family  of  algorithms  ineffective  for 
computing  combustion  problems  in  which  the  velocities  are  generally  low, 
but  where  the  flowfields  remain  strongly  compressible  because  of  heat 
release.  An  additional  low  speed,  compressible  flow  problem  that  is  of 
interest  to  the  present  authors  is  the  interaction  between  high  intensity 
radiation  fields,  including  both  high  power  laser  beams  and  focussed 
solar  radiation,  and  flowing  gases1 Our  purpose  in  the  present  paper 
is  to  identify  the  reasons  why  these  time-dependent  algorithms  fail  at 
low  speeds  and  to  devise  methods  for  enhancing  convergence  in  these 
regimes  so  that  they  can  be  used  effectively  for  these  additional 
applications . 

In  addition  to  extending  time-marching  methods  to  low  speed 
applications,  a  second  motive  for  studying  this  problem  is  to  enable  us 
to  understand  their  convergence  characteristics  more  thoroughly.  Besides 
increasing  the  range  of  convergence,  this  improved  understanding  may  also 
point  the  way  toward  improving  the  convergence  rate  in  those  regions 
where  the  methods  are  already  applicable.  In  the  present  paper,  we  study 
the  characteristics  of  time-marching  algorithms  in  this  light.  We 
address  the  primary  question  of  how  these  popular  algorithms  can  be 
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extended  to  low  Mach  number,  low  Reynolds  number  flows.  Throughout  our 
attention  is  limited  to  steady  flows. 

Methods  for  enhancing  low  Mach  number  convergence  in  inviscid  flows 
have  been  considered  by  several  previous  authors-^,  but  no  one  has  as 
yet  addressed  the  viscous  problem.  In  the  present  paper,  we  first  review 
this  previous  work  in  inviscid  flows  and  then  extend  these  procedures  to 
the  viscous  case.  Our  analysis  is  based  upon  implicit  algorithms  that 
use  central  differencing  in  space.  The  Mach-number/Reynolds-number 
regime  we  consider  ranges  from  high  subsonic  speeds  (Mach  numbers  of 
about  0.1)  down  to  incompressible  speeds  (Mach  numbers  around  10-6)  and 
from  inviscid  flow  (infinite  Reynolds  numbers)  down  to  highly  viscous 
conditions  (Reynolds  numbers  less  than  unity).  Appropriate  modifications 
to  the  implicit  algorithm  are  made  that  allow  a  single  unified  procedure 
to  give  efficient  convergence  over  this  entire  range.  Although  testing 
with  explicit  algorithms  has  not  been  attempted,  it  is  presumed  that 
these  procedures  will  also  provide  similar  improvements  in  convergence 
for  explicit  schemes  over  this  wide  regime.  The  philosophies  used  here 
should  also  prove  useful  for  developing  efficient  convergence  of 
flux-split  schemes  that  use  upwind  differencing. 

REVIEW  OF  CONVERGENCE  ENHANCEMENT  IN 
INVISCID  LOW  MACH  NUMBER  FLOWS 

The  convergence  rate  of  traditional  time-dependent  algorithms  for 
inviscid  flows  slows  down  as  the  Mach  number  is  reduced  because  of  the 
increasing  diversity  of  the  speeds  of  the  eigenvalues.  As  these 
eigenvalues  become  stiff,  both  explicit  and  implicit  algorithms  show 
slower  convergence  for  distinct,  but  related,  reasons.  Explicit  schemes 
slow  down  because  the  maximum  allowable  time-step  is  strictly  limited  by 
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stability  considerations  and  the  CFL  corresponding  to  the  slower 
eigenvalues  approaches  zero.  Implicit  approximately-factored  algorithms 
slow  down  because  the  factorization  introduces  an  optimum  CFL.  Slower 
convergence  is  observed  for  CFL's  above  or  below  this  optimum.  When  the 
Mach  number  is  low,  only  one  of  the  eigenvalues  (say,  for  example,  the  u 
or  u+c  eigenvalue)  can  be  kept  at  a  CFL  near  this  optimum  while  the 
others  are  far  from  the  optimum  and  convergence  slows  dramatically.  As 
is  shown  later,  fully  implicit  algorithms  that  use  direct  inversion  of 
the  complete  multidimensional  matrix  (no  approximate  factorization)  do 
not  show  a  convergence  slowdown  because  of  eigenvalue  stiffness.  They 
continue  to  show  rapid  convergence  rates  that  are  independent  of  Mach 
number.  The  reason  for  this  is  because  direct  inversion  methods  do  not 
exhibit  an  optimum  CFL,  but  continue  to  converge  more  rapidly  as  CFL  is 
increased. 

There  have  been  two  distinct  methods  proposed  for  circumventing  the 
convergence  slowdown  induced  by  eigenvalue  stiffness.  One  method  is  to 
use  time-derivative  preconditioning.  Early  studies^  0f  preconditioning 
showed  that  rapid  convergence  could  be  achieved  down  to  Mach  numbers  of 
about  0.01.  More  recently,  we  have  shown ^  that  preconditioning  allows 
Mach-number-independent  convergence  down  to  Mach  numbers  of  10~5t  but 
that  round-off  errors  begin  to  affect  the  maximum  convergence  level  at 
Mach  numbers  below  10"^.  In  this  regard  it  is  notable  that  it  is 
round-off  errors  in  the  pressure  (not  density  as  might  be  expected)  that 
eventually  prohibits  the  use  of  the  preconditioned  approach. 

The  second  method  for  circumventing  low  Mach  number  convergence 
difficulties  has  been  through  using  special  perturbed  forms  of  the 
equations  of  motion  that  are  valid  at  low  Mach  numbers.  Gustafsson^  has 


used  an  expansion  in  Mach  number  while  the  present  authors®  have  used  an 
expansion  in  Mach  number  squared.  Gustafsson’s  philosophy  was  to 
symmetrize  the  matrices  corresponding  to  the  Jacobians  of  the  flux 
vectors.  Our  Mach  number  squared  expansion  was  used  to  control  the 
disparity  in  the  eigenvalues  and  left  the  matrices  non-symmetric.  Our 
numerical  results  showed  effective  convergence  control  to  Mach  numbers  of 
10~®.  Calculations  at  lower  Mach  numbers  were  not  attempted,  but  we 
estimate  flows  to  10“®  or  ICM®  in  Mach  number  could  be  computed  before 
unacceptable  round-off  error  would  begin  to  decimate  the  results. 

Although  we  are  not  aware  of  similar  tests  of  Gustafsson' s  expansion 
procedure,  our  interpretation  of  his  approach  suggests  that  it  would  also 
be  effective  in  providing  Mach-number- independent  convergence  rates. 

Extension  of  these  low  Mach  number  methods  to  viscous  flows  has 
heretofore  not  been  attempted.  The  purpose  of  the  present  paper  is  to 
develop  a  method  that  allows  rapid  convergence  over  all  Reynolds  number 
regimes  without  sacrificing  convergence  rate  in  the  inviscid  case.  The 
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procedure  described  here  is  based  upon  the  perturbation  expansion  ! 

i 

procedure.  ; 

THE  CHOICE  OF  A  REPRESENTATIVE  PROBLEM  ; 

We  are  interested  in  assessing  the  convergence  rate  of  numerical  2 

H 

algorithms  over  a  broad  Reynolds-number/Mach-number  spectrum  ranging  from 

v 

v 

inviscid  to  highly  viscous  conditions  and  from  high  subsonic  speeds  to  ;> 

N 

very  low  velocities.  To  make  this  assessment  we  need  a  representative  ^ 

problem  that  will  possess  a  non-trivial  solution  over  the  entire  -i 

V 

Reynolds-Mach  number  regime.  As  such  a  problem  we  choose  the  flow  ^ 

> 

through  a  duct  with  volumetric  heat  addition  (Fig.  1).  The  heat  addition  J 

u 

causes  large  density  changes  such  that  even  in  the  presence  of  low  n 


velocities,  the  compressible  form  of  the  equations  must  be  used.  This 
problem  is  representati ve  of  typical  combustion  problems,  and  is  also 
analogous  to  the  radiation/gasdynamic  interactions^ >2  mentioned  above. 

In  the  absence  of  viscous  effects,  the  presence  of  a  spatially 
varying  volumetric  heat  source  prevents  the  solution  from  being  a 
trivial,  uniform  flow  and  instead  generates  an  inviscid  flowfield  that  is 
strongly  two-dimensional.  When  the  effects  of  viscosity  are  included, 
the  no-slip  condition  on  the  wall  ensures  that  the  velocity  profiles  will 
be  non-uniform.  At  very  low  Reynolds  number,  however,  fully  developed 
conditions  are  rapidly  established  and  the  introduction  of  a  heat  source 
again  provides  a  more  general  flowfield  at  either  low  or  high  Reynolds 
numbers.  For  viscous  flows  we  use  both  volumetric  heat  sources  and  heat 
addition  through  the  walls. 

CONVERGENCE  OF  THE  TRADITIONAL  TIME-DEPENDENT  PROCEDURE 
As  a  first  assessment  of  the  problem,  we  report  the  convergence  rate 
of  a  traditional  approximately-factored  scheme  with  no  provision  for  low 
Reynolds  number  or  low  Mach  number.  In  their  standard  form  the  equations 


4?  +  —  +  4-  *  4—  (v  +v)+-~(v  +  v  j 

dt  dx  dy  dx  v  xx  xy;  dy  v  yx  yy' 


where  the  vectors  Q,  E,  F,  Vxx,  VXy,  VyX  and  Vyy  are: 


pj2+p 


(e+p)u 


pv^+p 

(e+p)v 


3  =  partial  differential ;  p  =  rho 


Here,  p,  p,  u  and  v  are  the  pressure,  density  and  x  and  y  components  of 
the  velocity;  e  is  the  total  internal  energy  that  is  related  to  the  other 
parameters  by  e  -  p[e  +  'A(u2  +  v2)]  where  e  is  the  internal  energy.  The 
diffusion  terms  contain  the  viscosity,  p,  and  the  thermal  conductivity,  k. 
For  simplicity,  Stokes  hypothesis  has  been  used,  and  in  anticipation  of 
low  Mach  number  applications,  the  viscous  dissipation  has  been  dropped. 

The  speed  of  sound,  c,  is  given  as  c 2  -  yp/p  where  y  is  the  ratio  of 
specific  heat:. 

Traditional  approximately-factored  Euler  implicit  algorithms  applied 
to  Eqn.  1  lead  to  the  following1^ 


(I/At  +  ~-A  -  4-  R  4-) (I/At  +  4-B  -  4-  R  -~)AQ 

1  dx  dx  xx  dx1  v  dy  dy  yy  dy; 

-  -  [4^  +4^-4-  (V  +  V  )  -  4-  (V  +  V  )]  (4) 

Ldx  dy  dx  v  xx  xy'  dy  yx  yy'J  '  1 

where  superscripts  refer  to  the  time  step  level,  and  AQ  is  the  change  in 

Q  in  one  time  step  (AQ  -  Qn+^  -  Qn).  The  matrices  Rxx  and  Ryy  are 

appropriate  Jacobians  cf  Vxx  and  Vyy. 

The  rates  of  convergence  of  the  classical  approximate  factorization 
scheme  for  the  test  problem  given  in  Fig.  1  are  shown  in  Fig.  2  for  a 


u=  mu;  e  =  epsilon;  y  =  gamma;  a  =  delta 


wide  range  of  inlet  Mach  numbers.  At  high  subsonic  speeds  (0.7  Mach 
number)  the  scheme  converges  rapidly.  Convergence  to  machine  accuracy  is 
achieved  in  some  380  steps.  As  the  Mach  number  is  reduced,  the 
convergence  gets  much  slower  as  intimated  above.  At  a  Mach  number  of 
0.4,  the  convergence  has  already  slowed  down  by  more  than  a  factor  of 
two,  and  about  1000  iterations  are  needed  to  reach  machine  accuracy.  At 
Mach  0.1,  some  4500  iterations  are  required,  while  at  a  Mach  number  of 
0.01,  extrapolation  suggests  it  will  take  some  45,000  steps  to  reach 
machine  accuracy.  Because  of  the  difficulties  in  finding  optimum 
convergence  rates  at  such  slowly  converging  conditions,  the  results  shown 
here  were  all  computed  at  a  CFL  (based  on  u+c)  of  6.0.  At  the  lower 
Mach  numbers,  somewhat  faster  convergence  could  probably  have  been 
obtained  for  slightly  different  values  of  CFL,  but  the  rates  shown  here 
are  within  a  factor  of  two  of  their  optimum  value.  Clearly,  the  standard 
procedure  is  unacceptably  inefficient  at  low  speeds. 

The  physical  reason  for  these  convergence  difficulties  is  easily 
understood.  The  time-iterative  procedure  relies  upon  both  acoustic  waves 
and  particle  trajectories  to  propagate  errors  out  of  the  flowfield.  At 
low  Mach  numbers  the  acoustic  waves  make  many  trips  through  the  flowfield 
while  the  particles  are  traversing  it  a  single  time.  Approximate 
factorization  provides  most  rapid  convergence  when  the  individual  waves 
move  a  modest  fraction  of  the  distance  across  the  computational  domain  in 
a  single  time  step.  A  CFL  of  about  5  based  on  u+c  provides  this  optimum 
propagation  rate  for  the  acoustic  waves,  but  at  low  Mach  numbers  the 
correspond ing  CFL  based  on  u  is  so  small  that  the  particles  move  very 
small  distances  in  one  time  step  and  it  takes  many  steps  for  them  to 
traverse  the  flowfield.  For  this  reason  the  above  results  show  that  the 
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number  of  iterations  varies  approximately  inversely  with  the  Mach  number. 
Other  alternative  choices  of  CFL  do  not  provide  materially  faster 
convergence. 

Proof  that  it  is  the  approximate  factorization  errors  that  cause 
this  slowdown  in  convergence  is  easily  obtained  from  either  a  stability 
analysis  of  Eqn.  4,  or  from  its  direct  solution  without  approximate 
factorization.  If  we  define  the  amplification  matrix,  G,  as 

(G-I)q"  -  AQ  (5) 

the  Fourier  transform  of  Eqn.  4  gives  G  as  the  product  of  two  matrices, 

G  ■  K]~^K2  where. 


K,  -  I  +  i[S  +  S  “81  +  2[  ( 1-C  )4^  R  +  (1-C  )~?  R  1  +  T.c 
1  L  x  4x  y  Ay  J  L''  x'Ax^  xx  ;  y'Ay^  yyJ  AF 


K2  ■  1  -  SxSy  ESJ  >> 


(R  +  R  )  +  T.r 
'  *«  yx;  AF 


(6) 


Here,  Sx,  Cx,  Sy  and  Cy  represent  the  trigonometric  functions  of  the 
Fourier  modus  in  the  x  and  y  directions,  and  i  is  the  square  root  of 
minus  one.  The  term,  T/\p,  represents  the  errors  introduced  by 
approximate  factorizat  .on.  If  the  algorithm  is  solved  without 
factorization,  T^p  vanishes;  if  approximate  factorization  is  included,  it 
becomes , 

TftC  -  [iS  +  2 (1-C  R  ][ iS  +  2 ( 1-C  R  ]  (7) 

AF  L  x  Ax  v  x'Ax^  xx y  Ay  v  y'Ay6  yyJ  v  ' 


Although  these  expressions  for  the  ampl if ication  matrix  are  quite 
involved,  numerical  eigenvalues  are  easily  found  as  parametric  functions 
of  the  Fourier  wavenumbers  in  the  x  and  y  directions  for  specific  flow 
conditions.  Contour  plots  of  the  maximum  eigenvalue  of  G  obtained  from 
Eqn.  6  are  given  in  Fig.  3  for  flow  conditions  that  include  a  Mach  number 
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of  10"*  and  a  Reynolds  number  of  SO.  Figure  3a  shows  the  maximum 
amplification  factors  at  CFL  -  6  (based  on  u+c)  for  the  case  with 
approximate  factorization.  As  can  be  seen,  approximate  factorization  is 
stable,  but  the  eigenvalues  are  near  unity  everywhere. 

The  corresponding  stability  characteristics  without  approximate 
factorization  (T^p  -  0)  are  given  in  Fig.  3b.  Here,  we  have  used  a  CFL 
of  21,000  to  offset  the  low  Mach  number  and  the  amplification  rates  are 
much  less  than  unity  everywhere  suggesting  rapid  convergence.  Numerical 
experiments  with  the  fully  implicit  system  using  direct  inversion  rather 
than  approximate  factorization  verify  this.  In  general,  they  reach 
machine  accuracy  in  8  or  9  time  steps.  Direct-inversion ,  implicit 
procedures  eliminate  convergence  difficulties  at  low  Mach  numbers  and 
provide  a  very  robust  computational  algorithm.  Unfortunately,  the  CPU 
requirements  for  direct  solution  are  prohibitive  for  even  moderately 
refined  two-dimensional  grids  and  more  so  for  three-dimensional  problems, 
and  the  procedure  is  not  practical  for  routine  calculations. 

LOW  MACH  NUMBER  EXPANSION 

To  obtain  a  system  of  equations  that  is  valid  for  low  speed  viscous 
flows,  we  use  a  perturbation  expansion  similar  to  that  used  in  Ref.  8 
except  that  here  we  perturb  the  non-conservative  form  of  the  equations 
rather  than  the  conservative  form.  The  method  developed  in  Ref.  8  for 
inviscid  flows  is  unstable  at  moderate  and  low  Reynolds  numbers  and 
cannot  be  used  for  viscous  flows.  When  expressed  in  non-conservative 
form,  the  equations  are: 
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where  Q  -  (p,  u,  v,  p)1  and  the  Jacobians  A  and  B  are: 


The  vectors  V  ,  V  ,  V  and  V  are  as  given  in  Eqn.  3. 
xx  x.  yx  yy 

We  now  non-dimensionalize  these  equations  by  reference  values  of  the 
velocity,  density,  pressure,  and  temperature  which  will  be  denoted  by  a 
subscript  R,  and  by  a  reference  length  L.  Of  particular  interest  is  the 
momentum  equation,  and  here  we  consider  only  the  x-component.  In 
non-dimensional  form  it  becomes: 


-PR.  2  do 

Vr  P  dX 


:(V.T.) 


where  the  viscous  terms  (V.T.)  are  identical  in  form  to  those  in  Eqn.  8. 
The  ratio  of  reference  quantities  that  multiplies  the  pressure  gradient 
can  be  expressed  in  terms  of  a  reference  Mach  number, 


2m?  .  yM*  (11) 

PR  H 

where  y  is  the  ratio  of  specific  heats. 

Because  we  are  interested  in  flows  where  Mr  is  small,  we  define  a 
small  parameter,  e,  as, 

e  -  yMR2  (12) 

and  expand  all  parameters  in  a  power  series  in  e,  as. 
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(13) 


p  -  p0  +  epi  +  ... 

Inspection  of  the  non-dimensionalized  equations  shows  that  this 
substitution  affects  only  the  momentum  and  energy  equations.  The  1/e 
multiplying  the  pressure  gradient  in  the  momentum  equation  implies  the 
zeroeth  order  pressure  must  satisfy  the  relation: 

grad  p0  -  0  (14) 

This  indicates  p0  is  a  function  of  time  only.  For  steady  state  problems 
we  can,  with  complete  generality,  choose  p0  as  a  constant. 

Using  Eqn.  14  in  the  energy  equation  gives,  after  minor 
manipulation,  the  following  low  Mach  number  viscous  system: 

do  d  d  „ 
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where  V.T.  represents  the  appropriate  viscous  terms  as  given  in  Eqn.  3, 
and  Re  and  Pr  are  the  Reynolds  and  Prandtl  numbers  based  on  reference 
quantities.  For  simplicity,  we  have  dropped  the  zero  subscript  on  all 
variables  in  Eqn.  15  because  the  perturbation  expansion  causes  only  p]  to 
appear. 

Inspection  of  the  coupled  system,  Eqn.  15,  shows  that  there  are  five 
unknowns,  p,  u,  v,  T,  and  p-j ,  appearing.  The  density  and  temperature  are 
related  by  the  perfect  gas  law  which,  upon  expansion,  reduces  to  p0  -  pT. 
The  problem  with  Eqn.  15  is  that  the  time  derivatives  contain  relations 
for  both  p  and  T  along  with  u  and  v  but  there  is  no  provision  for 
updating  p]. 
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To  circumvent  this  problem,  we  replace  the  time  derivative  in  the 
continuity  equation  by  an  artificial  time  derivative  as  originally 
suggested  by  Chorin^  for  incompressible  flow.  Specifically,  we  replace 
the  continuity  equation  in  Eqn.  15  by, 

p0  dt  dx  dy  v  1 

where  6  is  a  constant  used  for  scaling.  With  this  form  of  the  continuity 
equation,  the  low  Mach  number  system  now  contains  p-j,  u,  v,  and  T  in  the 
time  derivative.  The  simplified  perfect  gas  relation  immediately  gives 
the  density.  Also,  we  note  that  the  final  three  equations  appear 
well-suited  for  viscous  effects  since  they  are  all  of  the  form: 

f£  +  C.T.  -  V2<|>  (17) 

where  <f>  represents  u,  v  and  T,  respectively,  and  C.T.  stands  for 
convective  terms. 

To  place  Eqn.  15  into  a  form  suitable  for  computation,  we  express  it 
in  conservation  form  using  Eqn.  16  as  the  continuity  equation.  Upon 
combining  these  we  obtain: 


V 


8  =  beta;  p  =  phi;  v  =  nabla;  r  =  gamma 


and 


Q*  -  (pi/8,  u,  v,  T) T 

E'  -  (pu,  pu2+p],  puv,  puT)T  (20) 

F*  -  (pv,  puv,  pv2+pi ,  pvT)T 

CONVERGENCE  OF  LOW  MACH  NUMBER 
SYSTEM  IN  INVISCID  FLOW 

To  test  the  perturbed  low  Mach  number  formulation  we  first  apply  it 
in  inviscid  flow.  Convergence  of  an  implicit  scheme  is  closely  related 
to  the  eigenvalues  of  the  flux  Jacobi ans.  For  this  low  Mach  number 
system,  the  eigenvalues  are: 

«  (u,  u,  (u+c')/2,  (u-c')/2)  (21) 

where  A 1  -  dE ' /dQ' ,  and 

c'2  -  u2  +  40T  (22) 

In  order  to  get  well-conditioned  eigenvalues,  we  choose  0  as  unity. 
Consequently,  the  flow  described  by  Eqn.  18  is  always  "subsonic",  an 
outcome  in  keeping  with  our  physical  expectations. 

Convergence  results  for  the  inviscid  model  problem  of  Fig.  1  based 
on  the  low  Mach  number  equations  are  given  in  Fig.  4  for  a  Mach  number  of 
10~3.  Identical  rates  of  convergence  are  obtained  at  Mach  numbers  of 
10~T  and  10"5.  The  expansion  has  completely  removed  the  dependence  on 
Mach  number.  The  convergence  observed  with  the  complete  equations  at  the 
high  subsonic  speed  is  observed  with  the  low  Mach  number  equations  at  all 
speeds. 


A  =  lambda 


CONVERGENCE  IN  LOW  MACH  NUMBER 
VISCOUS  FLOWS 

To  test  the  capability  of  the  present  low  Mach  number  expansion  in 
viscous  flows,  we  again  use  the  physical  problem  in  Fig.  1,  but  with 
no-slip  conditions  enforced  at  the  wall.  Representative  convergence 
rates  obtained  at  several  different  Reynolds  numbers  are  shown  in  Fig.  5. 
These  results  are  all  for  a  Mach  number  of  10“3.  At  high  Reynolds 
numbers  where  the  viscous  terms  are  small,  they  have  little  or  no  effect 
on  convergence.  As  shown  on  the  Figure,  the  rate  of  convergence  at  a 
Reynolds  number  of  1000  is  nearly  the  same  as  for  the  inviscid  case.  (In 
fact,  the  effect  of  viscosity  at  these  high  Reynolds  numbers  is  to  make 
the  rate  of  convergence  modestly  faster  than  in  the  inviscid  case.  This 
is  because  the  physical  diffusion  serves  to  dissipate  the  errors  in  the 
solution  slightly  more  rapidly.)  By  contrast,  the  rate  of  convergence  at 
low  Reynolds  numbers  is  much  slower  than  in  the  inviscid  case.  The  case 
at  a  Reynolds  number  of  10  in  Fig.  5  is  seen  to  converge  at  about 
one- fourth  the  rate  of  the  inviscid  case.  Flowfields  at  Reynolds  numbers 
less  than  10  failed  to  converge  entirely.  The  reasons  for  this  slowdown 
in  convergence  at  low  Reynolds  numbers  are  again  best  demonstrated  by 
investigating  the  stability  characteristics  of  the  low  Mach  number 
equations. 

Figure  6  shows  the  stability  characteristics  of  the  low  Mach  number 
for  the  three  Reynolds  number  cases  discussed  in  Fig.  5.  At  a  Reynolds 
number  of  1000,  the  stability  results  are  slightly  more  favorable  than 
for  the  inviscid  case  thereby  justifying  the  relative  convergence  rates 
for  the  two  conditions.  Values  over  much  of  the  wavenumber  domain  range 
between  0.90  and  0.95.  At  a  Reynolds  number  of  10,  the  amplification 
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factors  remain  stable,  but  are  very  near  unity  in  all  wavenumber  regimes. 
Typical  values  here  are  above  0.98.  These  larger  amplification  ratios 
explain  the  slow  convergence  that  was  observed  at  Re-10. 

The  reason  for  this  increase  in  amplification  factor  can  be 
understood  by  recalling  that  the  CFL  condition  of  inviscid  flow  is  joined 
by  a  second  dimensionless  quantity  in  viscous  flows  that  is  here  referred 
to  as  a  von  Neumann  number,  o  -  At/ (Re  Ax^).  The  slowdown  in  convergence 
at  low  Reynolds  numbers  is  because  the  von  Neumann  number  becomes  large 
and  the  approximate  factorization  error  in  the  diffusive  terms  drives  the 
ampl ification  factor  toward  unity  and  slows  convergence.  For  inviscid 
flow,  the  von  Neumann  numbers  are  zero.  At  a  Reynolds  number  of  1000, 
they  are  o  =  0.25,  while  at  Re  -  10,  o  has  increased  to  24.  Rapid 
convergence  requires  that  both  the  CFL  and  the  von  Neumann  numbers  be  of 
order  unity. 

The  most  obvious  way  to  control  both  the  Courant  and  the  von  Neumann 
numbers  simultaneously  is  by  changing  the  grid  size.  Figure  7  shows 
stability  results  for  Re=0.1  (which  is  a  worse  case  than  those  on  Fig.  6c) 
with  the  grid  spacing  being  increased  to  give  a  =  3.8  while  CFL  has 
been  fixed  at  8.  These  values  of  parameters  provide  amplification 
factors  that  are  considerably  farther  from  one,  and  promise  much  faster 
convergence  than  do  the  stability  characteristics  on  Fig.  6c. 
Unfortunately,  the  grid  size  required  to  attain  this  value  of  c  is  larger 
than  the  computational  domain  and  so  has  no  practical  utility.  The 
perturbation  expansion  does,  however,  give  us  an  alternative  method  for 
fixing  both  non-dimensional  parameters  that  is  just  as  effective,  but 
that  does  not  require  a  change  in  the  number  of  grid  points.  The  method 
for  achieving  this  is  described  in  the  next  section. 

a  =  sigma;  =  =  equals  approximately 


CONVERGENCE  CONTROL  AT  LOW  REYNOLDS  NUMBER 


For  the  viscous  flow  computations  shown  above,  the  scaling  parameter 
6  in  the  time  derivative  of  the  continuity  equation  was  set  to  unity  to 
provide  properly  scaled  eigenvalues  for  the  inviscid  terms.  As  the 
Reynolds  number  is  lowered,  the  inviscid  terms  become  less  important  and 
the  viscous  terms  begin  to  dominate.  At  these  low  Reynolds  number 
conditions,  it  is  no  longer  preferable  to  select  8  to  provide  properly 
scaled  eigenvalues  for  the  inviscid  terms.  Instead,  we  can  use  the 
parameter  8  to  enable  us  to  specify  both  the  von  Neumann  number  and  CFL. 

In  the  inviscid  limit  we  pick  the  step  size,  At,  by  an  appropriate 
CFL  number. 


CFL  - 


whereas  in  the  viscous  limit  the  appropriate  grouping  of  the  time  step  is 
the  von  Neumann  number. 


By  solving  Eqns.  23  and  24  for  At  and  equating,  we  find. 


2CFLAX 

(u+c‘) 
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Inserting  the  definition  of  c'  from  Eqn.  22,  and  solving  for  8  gives. 


n  u2  r,  2CFL  1 . 2  ■» i 

8  ”  4T  i^-R^u  'i 


where  Re^  1s  the  cell  Reynolds  number,  pRuRAx/pR.  Note  that  both  u 
and  T  are  non-dimensional  and  so  are  of  order  one. 


With  this  expression,  we  can  now  specify  both  the  von  Neumann  number 
and  the  Courant  number  and  then  use  Eqn.  26  to  compute  0.  Knowing  0  we 
can  then  compute  At  from  Eqn.  25  and  proceed  to  the  numerical  solution. 
Choosing  a-3.8,  CFl-8  and  computing  8  in  this  fashion  brings  the 
stability  results  of  Fig.  6c  to  values  shown  in  Fig.  8.  This  stability 
map  is  almost  identical  to  the  one  shown  on  Fig.  7  that  was  obtained  by 
changing  Ax,  and  again  promises  fast  convergence.  The  value  of  0 
required  to  obtain  the  results  on  Fig.  8  is  8  -  100. 

For  a  wide  range  of  Reynolds  numbers,  we  therefore  recommend  that 
the  parameter  0  be  chosen  according  to  the  relation: 

S  ■  mm  [1,  Wjjjjt.  -  1)  -!]}  (27) 

where  we  have  dropped  the  order  one  terms  from  Eqn.  26.  This  function  is 
plotted  as  a  function  of  the  cell  Reynolds  number  in  Fig.  9. 

Rates  of  convergence  for  several  values  of  Reynolds  number  using  the 
value  of  0  given  in  Fig.  9  are  shown  in  Fig.  10.  With  8  defined  as  in 
Eqn.  27,  the  rate  of  convergence  of  the  low  Mach  number  system  remains 
almost  unchanged  for  Reynolds  numbers  ranging  from  infinity  to  0.01. 
Clearly,  this  modification  makes  the  time- iterative  procedure  much  more 
robust. 

LOW  SPEED  HEAT  ADDITION  AT  VARIOUS  REYN0L0S  NUMBERS 
To  demonstrate  the  capability  of  the  low  Mach  number/low  Reynolds 
number  formulation,  we  present  flowfield  results  for  the  heat  addition 
problem  described  in  Fig.  1.  A  series  of  plots  showing  tempe’n*'"'0  and 
velocity  contours  for  a  fixed  volumetric  heat  addition  and  an  inlet  Mach 
number  of  1  x  10-5  are  given  in  Fig.  11  for  Reynolds  numbers  of  infinity 


(inviscid),  5000,  500,  50,  and  0.05.  The  qualitative  effects  of 
viscosity  are  easily  recognized  in  these  plots.  For  inviscid  flow  the 
temperature  contours  are  convected  out  of  the  flowfield  with  no 
diffusion,  and  the  slip  boundary  condition  shows  no  boundary  layer  on  the 
wall  in  the  velocity  contour  plot.  At  a  Reynolds  number  of  5000,  a 
modest  diffusion  of  the  temperature  from  the  hot  to  the  cold  fluid  is 
noted,  and  a  thin  boundary  layer  exists  near  the  wall.  Reducing  the 
Reynolds  number  to  500  and  then  to  50  shows  the  increasing  effects  of 
heat  diffusion  and  a  thicker  boundary  layer  as  the  effects  of  diffusivity 
spread  over  nearly  the  entire  flowfield.  Jumping  to  a  Reynolds  number  of 
0.05,  the  effects  of  convection  are  now  completely  absent  in  the 
temperature  contours.  The  temperature  contours  give  an  undistorted  image 
of  the  volumetric  heat  source.  The  velocity  field  at  this  low  Reynolds 
number  condition  becomes  nearly  fully  developed  after  the  peak  in  the 
heat  source  is  passed  and  shows  no  reminiscence  of  boundary  layer  flow. 
The  value  of  B  used  for  each  case  is  given  in  the  Figure. 

Additional  flowfield  solutions  showing  the  effect  of  heat  addition 
through  the  walls  are  given  in  Fig.  12.  Here,  the  temperature  of  the 
wall  is  linearly  increased  from  T-l  at  the  inlet  to  T-6.67  at  the  exit, 
and  the  flow  is  heated  through  the  boundary  layer.  In  the  inviscid  case 
the  flow  is  not  affected  by  wall  temperature  and  only  viscous  results  are 
shown.  At  a  Reynolds  number  of  5000,  both  temperature  and  velocity 
contours  show  thin  boundary  layers  near  the  wall.  The  temperature  in  the 
middle  of  the  duct  is  unaffected  by  heating.  At  a  Reynolds  number  of  50, 
the  boundary  layers  are  much  thicker  and  the  heat  has  diffused  to  the 
middle  of  the  duct.  At  a  Reynolds  number  of  0.05,  the  diffusion  rates 
are  high  enough  relative  to  the  convective  speed  that  the  temperature  is 


uniform  in  y.  The  monotonic  increase  in  temperature,  however,  causes  the 
flow  on  the  centerline  to  continue  to  accelerate. 

This  spectrum  of  flow  conditions  ranging  from  inviscid  to 
viscous-dominated  flow  shows  that  the  modified  equations  provide 
qualitatively  correct  variations  in  the  velocity  and  temperature  fields 
as  the  magnitude  of  the  Reynolds  number  is  varied.  The  fact  that  all 
calculations  were  obtained  with  the  same  algorithm  and  that  all  showed 
nearly  identical  rates  of  convergence  is  an  indication  of  the  broad 
capabilities  of  the  approximately  factored,  implicit  time-dependent 
algorithm  when  proper  mathematics  and  physics  are  incorporated. 

DISCUSSION  AND  CONCLUSIONS 

Time-marching  algorithms  are  widely  used  for  the  solution  of 
transonic  and  supersonic  flows  and  for  problems  containing  embedded  shock 
waves.  These  methods  are  equally  applicable  to  inviscid  and 
high-Reynolds-number  viscous  calculations.  The  algorithms  provide 
accurate,  efficient  solutions  in  inviscid  flows  because  they  treat  the 
convective  terms  appropriately.  The  addition  of  weak  (high  Reynolds 
number)  viscous  effects  to  this  properly  formulated  inviscid  procedure 
has  little  or  no  effect  on  convergence  because  the  convective  terms  still 
dominate  the  flow.  The  convergence  rate  of  the  algorithm,  however, 
becomes  notoriously  slow  in  regions  where  the  Mach  number  is  low. 
Experience  at  very  low  Reynolds  number  conditions  is  very  limited 
(largely  because  of  problems  with  the  convective  terms),  but  all  evidence 
suggests  a  similar  slowdown  in  convergence  at  low  Reynolds  numbers. 

These  limitations  can  seriously  impair  convergence  rates  in  high 
speed  flows  because  of  local  low  speed  regions  adjacent  to  stagnation 
points  or  because  of  near  wall,  low  Reynolds  number  effects  in  boundary 
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layers.  Further,  the  limitation  to  high  speed  flows  precludes 
computations  of  important  new  low  speed  problems  such  as  combustion 
applications  in  which  the  strong  heat  addition  requires  the  use  of  a 
compressible  formulation  even  though  the  Mach  numbers  would  suggest 
incompressibility.  The  purpose  of  the  present  paper  has  been  to  identify 
and  remove  these  limitations  on  convergence  so  the  time-march ing 
procedure  can  be  applied  over  a  much  broader  regime. 

Although  previous  papers  have  studied  the  inviscid  low  Mach  number 
problem,  none  have  considered  the  viscous  problem.  Simple  computational 
experiments  show  that  neither  previous  preconditioning  procedures  nor 
previous  perturbation  expansion  procedures  provide  approaches  that  are 
effective  in  viscous  flows.  Accordingly,  a  new  low  Mach  number  expansion 
procedure  was  developed  from  the  non-conservative  equations  that  is 
applicable  to  both  viscous  and  inviscid  conditions.  This  expansion 
procedure  provides  a  mechanism  for  properly  scaling  the  eigenvalues  of 
the  convective  terms  thus  ensuring  fast  convergence  in  inviscid  flows  at 
all  Mach  numbers.  It  also  provides  a  mechanism  for  keeping  both  the  von 
Neumann  number  and  the  Courant  number  of  order  one  at  very  viscous 
conditions,  thereby  ensuring  rapid  convergence  at  low  Reynolds  numbers. 

The  resulting  domain  of  convergence  of  the  low  Mach  number,  low 
Reynolds  number  procedure  is  summarized  on  Fig.  13.  Rapid,  effective 
convergence  of  the  procedure  is  observed  for  Reynolds  numbers  ranging 
from  infinity  down  to  0.05  and  from  Mach  numbers  between  0.3  and  10“®. 
Rapid  convergence  is  here  defined  as  convergence  rates  within  20-30%  of 
that  obtained  in  inviscid  flows  at  Mach  0.7.  As  indicated  on  the  figure, 
the  modified  procedure  gives  slower  convergence  over  a  small  corner  of 
the  Reynolds  number/Mach  number  domain.  For  comparison,  the  much  smaller 
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region  of  convergence  of  the  original  time-marching  procedure  for 
compressible  flow  is  shown  in  the  upper  right-hand  corner  of  the  Re-Mach 
number  domain.  Clearly,  the  modified  procedure  based  on  the  perturbation 
expansion  makes  the  classical  compressible  algorithms  much  more  robust 
over  a  much  wider  range  of  variables.  In  particular,  this  low  Mach 
number  procedure  represents  an  effective  numerical  algorithm  for  the 
computation  of  low  speed  flows  with  strong  heat  addition. 
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